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Abstract 

Here we review a hybrid lattice Boltzmann algorithm to solve the equations 
of motion of cholesteric liquid crystals. The method consists in coupling a 
lattice Boltzmann solver for the Navier-Stokes equation to a finite difference 
method to solve the dynamical equations governing the evolution of the liq- 
uid crystalline order parameter. We apply this method to study the growth 
of cholesteric blue phase domains, within a cholesteric phase. We focus on 
the growth of blue phase II and on a thin slab geometry in which the domain 
wall is flat. Our results show that, depending on the chirality, the growing 
blue phase is either BPII with no or few defects, or another structure with 
hexagonal ordering. We hope that our simulations will spur further exper- 
imental investigations on quenches in micron-size blue phase samples. The 
computational size that our hybrid lattice Boltzmann scheme can handle 
suggest that large scale simulations of new generation of blue phase liquid 
crystal device are within reach. 
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1. Introduction 

The lattice Boltzmann (LB) algorithm [U [2] is a powerful method to solve 
the Navier-Stokes equations of ideal and complex fluids, which, due to its 
conceptual simplicity and ease to code, provides an attractive alternative to 
other methods more commonly adopted in computational fluids dynamics, 
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such as e.g. finite elements algorithms. In the last couple of decades, the 
LB method has evolved to a powerful tool to study a variety of problems 
in complex fiuids [3ll2|5l[6l[71El|9l[l0l[IIl[l2lini[Il]. In this respect, 
particularly successful avenues have been, among others, the study of binary 
fiuids and of liquid crystals. 

The original LB method to solve the hydrodynamic equations of mo- 
tion of complex fiuids typically consisted of introducing extra distribution 
function, which evolved via an LB dynamics determined by an appropri- 
ate collision kernel. These distribution functions control the dynamics of 
conserved or non-converved order parameters which is coupled to the fiuid 
momentum, itself governed by the primary LB distribution function. The 
Chapman-Enskog expansion of such dynamics then gave back the contin- 
uum hydrodynamic equations for the coupled system. However, on top of 
introducing some systematic, albeit small, errors, this "full LB" approach 
had the drawback of requiring a large memory to store all the distribution 
functions for the whole lattice. This turned out to be particularly severe for 
the case of liquid crystal, where the order parameter (in the Beris-Edwards 
model [15]) is a traceless and symmetric tensor, and therefore required the 
introduction of five extra sets of distribution functions. This leads to high 
memory requirements which ultimately limits the size of the lattice to be 
studied. The new generation of LB studies for complex fiuids has a growing 
need to be able to cope with larger and larger systems, to make best use of 
the potentially very good scalability of parallel LB codes. As a result, new 
hybrid algorithms have been coded and deployed, for both binary fiuids and 
liquid crystals [H HH [16], where the LB algorithm is dedicated to solve the 
(forced) Navier-Stokes equation for momentum alone, and is coupled to a 
standard efficient finite-difference solver (e.g. a predictor corrector) for the 
order parameter dynamics. Hybrid LB simulations of binary fiuids [9] and 
of active gel and active liquid crystal hydrodynamics and rheology pT] [12] 
have validated the hybrid algorithm through a comparison with full LB sim- 
ulations, and have shown that this framework is potentially very fiexible and 
robust. Here, after briefiy reviewing the equations of motion of liquid crystal 
hydrodynamics and the hybrid LB approach which we use to solve them, 
we apply our algorithm to study the growth of domains in cholesteric blue 
phases (BPs). 

BPs are a spectacular example of a soft solid, made up by a spontaneously 
occurring network of disclination lines in a cholesteric liquid crystal [T7]. BPs 
appear very close to the transition between the isotropic and the cholesteric 
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phase. In a liquid crystal in the cholesteric phase, the director field, which 
quantifies the local direction of molecular ordering in the fluid, has the ten- 
dency to twist around an axis, which is that of the cholesteric helix. Close 
to the transition to the isotropic phase, a simple helical order is frustrated 
and it is more advantageous for the director fleld to rotate in a helical fash- 
ion around any axis perpendicular to a line - this director field pattern was 
named a "double twist cylinder". Mathematically, it is impossible to patch 
up different double twist cylinders without creating defects in between, and 
this is what gives rise to the disclination network making up BPs. 

Without an electric field, three different kinds of BPs are found. BPI and 
BPII are cubic phases characterised by a regular disclination line lattice: in 
BPI the double twist cylinders are arranged in a simple cubic lattice, whereas 
in BPII they form a body-centred cubic lattice. BPIII, or the "blue fog", 
is not a cubic phase, and its structure is not fully understood to date: the 
dominating view is that it is a structure made up by double twisted cylinders 
arranged in an irregular way [T7]. 

As the period of the BP network and the timescales of response to an ap- 
plied perturbation are in the sub-/xm and sub- /is range respectively, these soft 
materials are promising candidates for tunable photonic crystals and a new 
generation of fast liquid crystal devices. Recent experiments [18] managed 
to stabilize BPs over a strikingly wide temperature range of 50 K, compared 
to the preceedingly narrow temperature range of about 1 K, putting these 
technological advances now within reach. However, in order for this excit- 
ing potential to be fully realised, our understanding of BPs needs to become 
as quantitative as the one we have for conventional nematic liquid crystal 
displays. 

From a more fundamental physics points of view, BPs are remarkable 
materials. A series of very interesting experiments EHl EI] have shown 
for instance that BP droplets facet when they nucleate and grow inside an- 
other isotropic fluid. BP droplets, like other fluids, may also wet a surface, 
but due to their elastic properties their surface shows steps and may re- 
construct, like that of a solid. According to the current understanding, the 
underlying disclination meshwork is primarily cause of this highly non-trivial 
phenomenology, as it behaves like an elastic network. 

Existing simulations of BPs [22| [23| IMf [25] have signiflcantly extended 
our quantitative understanding of their physics, which flrst rested on semi- 
analytical approximations (see e.g. [2S])- For instance, obtaining the shape 
of the phase diagram, with BPI and BPII appearing in the experimentally 
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observed order upon increasing the chirality [221 E3] , as well as understanding 
the presence of anomalous electrostriction in BPI under an electric field [25] , 
were only possible with extensive simulations of these structures. However, 
these simulations have thus far been limited to one unit cell, within which 
several disclination cores are present and require a fine enough discretisation 
to be correctly resolved. By necessity, this approach leaves out a number of 
physically relevant supra-unit cell phenomena such as the possibility of large 
scale lattice reconstructions (perhaps induced by shear or by an applied field), 
and the appearance of defects or spatial inhomogeneities in the BP lattice. 
Furthermore, the length scale covered by simulations is much smaller than 
the ones relevant for either the domain growth experiments or BP devices. As 
a first step to model more realistic situations, we present here supra-unit cell 
large scale 3D simulations of the domain growth of BP domains in cholesteric 
and isotropic liquid crystals. 

Our hybrid lattice Boltzmann results suggest that the physics of BP do- 
main growth is highly non-trivial. Focussing on a thin slab geometry, we 
observe below qualitatively different growth dynamics of BP domains inside 
cholesterics, for different values of the chirality. For large values of this pa- 
rameter, we find that BPII domains growing inside a cholesteric evolve into 
a different disclination network, which possesses local hexagonal symmetry. 
Further simulations in an isotropic fluid show the same transition, although 
the growth kinetics of the new BP phase is different. To aid comparison 
with potential future experiments, we also visualize the disclination network 
and simulate the appearance of the sample under polarized light. The fea- 
sibility of such large scale simulations is promising in view of further future 
applications, e.g. to device modelling. 

2. Model and methods 

2.1. Equations of motion 

Following the Beris-Edwards model for liquid crystal hydrodynamics [13], we 
describe the BPs by a traceless, symmetric, second rank order parameter 
tensor, Q. The equilibrium thermodynamic properties are determined by a 
Landau - de Gennes free energy JF, whose density / is, 

•f ~ " 3)Qal3 ^QapQfB-yQ-ya -\ f^iQap)'^ 

+ y (CaiJ^lQiS/? + 2goQa/3) • (1) 
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(Notice that in our notation Greek indices denote Cartesian components and 
summation over repeated indices is implied.) In Eq. [l]y4o is a constant, K is 
an elastic constant (we adopt the one-constant-approximation), go is 27r/p, 
where p is the pitch of the cholesteric helix, while 7 is a control parameter 
linked to the temperature for thermotropic liquid crystals, or concentration 
for lyotropic ones. Increasing 7 leads to an increase in the average magnitude 
of order. The local magnitude of order, q{r), is the largest eigenvalue of 
Q. This quantity is also used to identify disclination lines: when the local 
magnitude of order falls below a predetermined threshold, we identify that 
lattice point as constituting part of a disclination. This prescription is very 
easy to implement and allows an accurate determination of the disclinations. 
Changing the numerical value of the threshold generally leads to a change 
in the thickness of the disclination tubes. For our study, we have typically 
chosen a threshold of g = 0.19 for defect/disclination identification. The 
Beris-Edwards equations for the evolution of the Q-tensor, which we aim to 
solve, are [U] 

AQ = r(^ + |Tr(g)i)=rH. (2) 

Here F is a collective rotational diffusion constant, while ^ indicates the 
functional derivative with respect to the tensor order parameter. Tr stands 
for trace, and Dt is the material derivative for rod-like molecules [15] : 

AQ = (9t + M- V)Q-S(W,Q) (3) 
S(W,Q) = (eD + u;)(Q + I/3) (4) 
+ (Q + I/3)(eD-^)-2^(Q + I/3)Tr(QW), 

where D and uj are the symmetric and the anti-symmetric part respectively 
of the velocity gradient tensor Wap = dpUa [iSl HO] , u being the velocity field. 

Eq. |2] may be mathematically derived from an underlying microscopic 
model via a Poisson brackets approach [15]. Physically, this equation means 
that the system tries to evolve, in the absence of any flow or backflow (^2 = 0), 
so as to minimise its free energy. This is ensured by the presence of the molec- 
ular field, H. The presence of a non-zero u requires one to substitute the 
usual partial derivative dt with the material derivative Dt-, which describes 
advection by the fluid velocity (the term u ■ VQ), and also includes a fur- 
ther coupling S(W, Q) between the velocity gradient tensor and the order 
parameter, which arises due to the tensorial nature of the latter. Physically, 
S(W, Q) appears because the order parameter distribution can be both ro- 
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tated and stretched by flow gradients [15]. Mathematically, S(W, Q) con- 
tains a (phenomenological) mixture of upper and lower convected derivatives. 
This is similar to what is done in other hydrodynamic equations for poly- 
meric fluids, such as the Johnson- Segalman model [27]. The term ^, which 
may be also viewed as a parameter controlling the weights of the lower and 
upper convected derivatives, is related to the aspect ratio of the molecules. 
This parameter controls whether the director field is flow aligning in shear 
flow > 0.6), creating a stable response, as opposed to flow tumbling, which 
gives an unsteady response < 0.6). 

The fluid velocity, m, obeys the continuity equation of an incompressible 
fluid, and the Navier-Stokes equation, 

p{dt + updp)uc, = dfsiHaf}) + r]df3{daUfi + dpUa (5) 

where p is the fluid density, r] is an isotropic viscosity. The stress tensor IIq,^ 
consists of a symmetric part cTq,^ and an antisymmetric contribution Tq,^, 

^afS = (^af3 + TafS (6) 
^0/3 = —PoSal3 + '^^{Qal3 + lSaf3)Q^eH^e (7) 

df 



la/3 



= Qa^H^p — HayQ"ff3, (8) 



where Pq = pT — / is an isotropic pressure [15] and Q-^/u,i3 = dpQ.^ 

In our simulations, the coupling to hydrodynamics via a non-trivial pres- 
sure tensor can be switched off, essentially by imposing a constant zero ve- 
locity proflle. In this way the effects of flow and backflow may be unambigu- 
ously pinpointed. In our case, backflow does not dramatically modify the 
kinetic pathway through which BP domains grow, but it renders the dynam- 
ics faster. Similar effects have been seen on the switching dynamics and on 
the rheological properties of BPs [211 [25] . 

The equilibrium phase diagrams of BPs are commonly expressed as a 
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function of the chirality k and reduced temperature r, defined [221 ESI [2S] as 



K 
T 

We selected r = and varied k in our simulations (reported in the Results 
section). For the other kinetic parameters, we chose F = 0.3, ^ = 0.7 (hence 
a flow aligning BP) and r] = 5/3 (these choices are suggested by previous 
experience in blue phase numerics 125] ) . Time and space is measured in 
simulation units in what follows. 

2.2. The hybrid lattice Boltzmann algorithm 

In order to solve Eqs. |2]and[5| as mentioned in the Introduction, we use here a 
hybrid lattice Boltzmann algorithm. The idea of the algorithm is simple. We 
observe that the coupling between the velocity field and the order parameter 
equation (see Eq. |2]) is via the material derivative term, which requires 
both the velocity and the velocity gradient fields. On the other hand, the 
order parameter field affects the Navier-Stokes equation through the pressure 
tensor 11^^. Our hybrid lattice Boltzmann approach consists in solving Eq. [2] 
via a finite-difference predictor-corrector algorithm, while the LB algorithm is 
devoted to the integration of Eq. |5| The order parameter and velocity fields 
are (sequentially) updated at every time step via these algorithms. The LB 
step requires as an input the order parameter field (hence the pressure tensor 
and its divergence), which is provided by the finite difference solver. On the 
other hand, the LB algorithm updates the velocity field which is then in turn 
required by the finite difference scheme to further evolve the dynamics of the 
order parameter. Due to the limited required coupling between LB and finite 
difference algorithms, this hybrid approach could be generalised to a variety 
of hydrodynamic equations of motion involving the dynamics of an order 
parameter besides the Navier-Stokes equations [HI ttH US]- As mentioned 
above, a similar hybrid approach has indeed already been successfully used 
for a binary fluid as well |9i il6j . 

With respect to a full LB approach [281 HOj, the primary advantage of 
this hybrid method is that it will allow simulations of larger systems as it 
involves consistently smaller memory requirements. Indeed, while in a full 
LB treatment one has to store 6 sets of 15 distribution functions at any 
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(10) 
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lattice point (if we choose the 3DQ15 velocity vector lattice pQ as we do 
here), in the hybrid approach just one set of distribution functions plus the 
five independent components of the Q tensor are needed. The hybrid method 
should also be numerically more stable since it avoids the error term arising 
in the Chapman-Enskog expansion used to connect the LB model to the 
order parameter evolution equation in the continuum limit [28J. Finally, it 
should be possible to include noise in a conceptually simpler way as one 
can in principle couple a noisy LB dynamics, for which precise prescriptions 
exist [29], for the Navier-Stokes equation with a Langevin dynamics for the 
order parameter(s) evolution. In the present work, however, we neglect noise. 
Thus our treatment couples a mean-field dynamics of the cholesteric order 
to a noiseless classical fluid. 

Let us describe first how to integrate the Navier-Stokes equations with a 
slightly modified standard Lattice Boltzmann algorithm ^28j. This is defined 
in terms of a single set of partial distribution functions, the scalars fi{x), 
that sum on each lattice site x to give the density. (Note that / without 
Latin index is the free energy density from Eq. Each fi is associated 
with a lattice vector Cj [281 ITO] . We choose a 15- velocity model on the cubic 
lattice with lattice vectors: 

(0,0,0) (11) 
(±1,0,0),(0,±1,0),(0,0,±1) (12) 
(±1,±1,±1). (13) 

The indices, i, are ordered so that i = corresponds to ef\ i = I,-- - ,6 
correspond to the e^^^ set and i = 7, - ■ ■ , 14 to the e^^^ set. 

Physical variables are defined as moments of the distribution functions: 

P=^fu pUa = ^ fiCia- (14) 

i i 

The distribution functions evolve in a time step At according to 

Uix+eAt, t+At)-fi{x, t) = Y P/^^^' ^' ^^*}) + ^f'^^ + ^ + ■ 

(15) 

This represents free streaming with velocity Cj followed by a collision step 
which allows the distributions to relax towards equilibrium. The /*'s are 
first order approximations to fi{x + eidt,t + dt), and they are obtained by 



= 

-(1) 
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using At Cfi{x, t, {fi}) on the right hand side of Eq. ( 15 ). Discretizing in this 
way, which is similar to a predictor-corrector scheme, has the advantages that 
lattice viscosity terms are eliminated to second order and that the stability 
of the scheme is improved [28] . 

The collision operators are taken to have the form of a single relaxation 
time Boltzmann equation, together with a forcing term 

C^x, t, {/,}) = t) - (f , t, {/,})) + p,{x, t, {/,}), (16) 

The form of the equations of motion follow from the choice of the moments 
of the equihbrium distributions f^'^ and the driving terms pi. Note that /^^^ 
is constrained by 

5Z ^i'^ ^ P' 5Z •^i"^^*" ^ X] = -cr^/j + pUo^Up, (17) 

i i i 

ff^eiaCipei^ = ^ {UaSf3^ + Uf36aj + U^^ap) (18) 

i 

where the zeroth and first moments are chosen to impose conservation of 
mass and momentum. The second moment of /^^^ is determined by a^p, 
which is the symmetric part of the stress tensor Hap, and does not include 
either the double gradient term, daQ-yu , or the fSap contribution from 
the isotropic pressure. The constraint on the third moment is necessary to 
get an isotropic Navier-Stokes equation via the Chapman-Enskog expansion 
of Eq. |15| (see e.g. [28j). 

The divergences of and of daQ-yu an^ instead enter effectively as 
a body force and constrain the first moment of the driving terms pi: this 
ensures that spurious velocities are eliminated, or greatly reduced, as we 
show below. The constraints on the pj's are [30] : 

y^Pi = 0, y^^PiCia = dpTaH - dp (daQju ] + ^af = (19) 

i i V (yQ'ru,i3j 

y^PieiaCi/j = 0, y^^Pieigei/Bej-y = 0. (20) 



This scheme was inspired by the one used in Ref. [31] to reduce spurious 
velocities in a full LB binary fluid simulation. In the case of liquid crystals 
spurious velocities are equally eliminated by introducing the divergence of 
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Qap = —daQ-yuQ^ 1- f^ap- To 866 this, W6 transform dpQap in th6 folfowing 

way: 



df 



- da 



+ 



Of 



dQ 



daQ 



+ daf 



df 



-daQ 



df ^ df 



dQyu dQ 



dQfu,/3 
df^Q^i/ H^pdf^Q^p^ 



(21) 



wher6 W6 introduc6d H^y in th6 last lin6 b6caus6 th6 Q t6nsor is trac6less and 
th6r6for6 Q^^ = 0. W6 S66 that if W6 dir6ctly ins6rt Eq. 21 as a body force 



this term vanishes in equihbrium as it is proportional to the molecular field. 
In this way spurious velocities, i.e. non-zero velocities in equilibrium due to 
the different discretisation of the pressure tensor and the molecular field, are 
in principle avoided. We have verified that this is the case numerically. 



Conditions Eqs. 17 - 20 are satisfied by writing the equilibrium distribu- 
tion functions and forcing terms as polynomial expansions in the velocity, 

(22) 

where s G {0,1,2} ■<=^ ^ G {0,1,3} identifies separate coefficients for 
different square absolute values of the velocities. 

The coefficients in the expansion are determined by the requirements that 
all the constraints enumerated above, Eqs. T7|l8||19 and 20, are fulfilled. A 
possible choice, which we adopt, is given below: 



A 



= pL 

' 10' 
B2 = pI2A, 

_P_ 
24' 
P_ 
16' 



C2 = 
D2 = 

E2al3 

P2 = 



Ai = A2, 
Bi = 8B2, 
Ci = 2C2, 

Di = 8D2 



1 

24' 



1 

fe ^"""^ 3 
Pi = 8P2. 



= p - 14^2, 



El 



a/3 



2p 
3 ' 



8E2a(i 



(23) 
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Clearly the stress tensor Ua/B is a function of Qaf3 and, in the hybrid approach, 
the coefficients E2a/3 and Eia/s of the /f^'s are computed by using the solution 

(via finite difference methods) of the coupled Eq. ([2]). This differs from the 
fully LB treatment of liquid crystals [28| ITO] . 

2. 3. Rendering of the optical pattern 

Due to their anisotropic structure, liquid crystal molecules are optically ac- 
tive and cause polarisation rotations and phase-shifts in the electric field 
components of transmitted polarised light. A common experimental tech- 
nique is to observe the transmission pattern under a microscope using a 
crossed-polariser geometry. Every molecule along the path of the incident 
beam acts as a retarder. Thus the transmission pattern can be regarded 
as an overall effect of the constituting molecules. It depends on the local 
orientation of the molecules, i.e. on the director field, the refractive-index 
anisotropy and on the shape of the domain. 

While it is in general not trivial to infer the local director field from the 
transmittion pattern, it is rather simple to simulate the latter from a given 
director field. To this end, we have employed the Miiller matrix technique [321 
1^ . which simulates the light transmission signal observed in a micron-size 
sample under a pair of crossed polarisers. We calculated the polarised optical 
texture corresponding to the instantaneous director field of the growing BP 
domain and its surrounding environment. In our approach the director field 
d{r) is defined as the normalized eigenvector related to the largest eigenvalue 
of the order parameter tensor Q. According to Stokes the polarisation state of 
the light can be conveniently described by a set of four parameters, combined 
in the 4-component Stokes vector S = {Sq, Si, S2, S3). Its first component 5*0 
is proportional to the intensity of the light. 

In order to simulate the retardation of a liquid crystal droplet one assumes 
the droplet of size L to consist of N equally thick liquid crystal layers with 
thickness h = L/N. The effect of each layer situated at the site r onto the 
Stokes-vector 5* is then described by the Miiller matrix 



/ 1 










(24) 



where 



Cb = cos2P{r), Sb = sin2/9(r), Cd = cos 6{r), Sd = sin 6{r) 



(25) 
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and 



a{r) = arccos ((/^(r)) , (26) 
(3{r) = arctan(4(f)/o?j,(f)), (27) 
2TTq{r)h ( UoU^ \ 

Vv^o^i^ Q;(r\) + rig cos^ a (rj J 

The quantity is the local scalar order parameter defined as the largest 
eigenvahic of the order parameter tensor Q, whereas A gives the wavelength 
of the incident light. We set the ordinary and extraordinary refractive indices 
to Uq = 1.5 and rig = 2.0 respectively. 

The above definitions of the angles a and /3 are correct for a light beam 
propagating along the z-direction. We assumed this throughout our anal- 
ysis, but the adaption to other situations is straightforward. Generally, a 
is the angle between the direction of the light beam and the local director 
field, while (3 measures the angle between the projection of the local director 
field onto the coordinate plane perpendicular to the beam direction and a 
coordinate axis in that plane, in our case the x-axis. The crossed-polariser 
geometry is realized by two different Miiller matrices of the type 



/ 1 COS0 sin0 \ 

cos cos^ (f) cos sin 

sin (f) sin cos (j) sin^ 

\ 0/ 



(29) 



The parameter defines the angle between the polariser or analyzer and 
a coordinate axis perpendicular to the beam direction. For a right angled 
crossed-polariser setup one assumes for instance 0j„ — for the polariser and 
(pout = 7r/2 for the analyser. 

The total effect of the liquid crystal droplet on a Stokes vector Sin can 
be formally expressed by a matrix product of consecutive Miiller matrices, 
following the path of the hght beam: 

Saut{x, y) = Pout n^=o M{i = x/hj = y/h, k) PinSin{x, y). (30) 

The matrices M {i.j, k) shall be understood as the discretised version of M{r). 
As we are interested in the general appearance of the sample under poly- 
chromatic polarised light, we performed this operation for slightly different 
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wavelengths in the magnitude of the unit cell size of the BPs, which was 16 in 
lattice units, and weighted the results accordingly [3l]. We simulated three 
different components with wavelengths A = 16, 18 and 20 in lattice units and 
assigned them the weights w = 0.2,0.6 and 0.2 respectively. This procedure 
can only be regarded as a simple approach to model white light. Nevertheless 
it proves to be sufficient as the results for different components with similar 
wavelengths differ only slightly in their brightness and only insignificantly in 
their general transmission patterns. The results for the first component 5*0 
of the Stokes vector S are displayed in grayscale in Fig. [ij 

3. Results 

We now discuss the results obtained by applying our hybrid lattice Boltzmann 
algorithm to study the growth of blue phase domains inside cholesteric liquid 
crystals. To better analyse the growth dynamics, we restricted ourselves to 
the case of a thin slab of liquid crystal, and assumed periodic boundary 
conditions along the small direction, which we took along the z axis. This 
corresponds to assuming that the domain wall is locally straight and should 
approximate well the experimental situation in which the growing blue phase 
droplets are very large with respect to the unit cell size, so that locally they 
can be viewed as basically planar. In this geometry, which only requires few 
(16-32) lattice points along the z direction, we can study the time evolution of 
the order parameter and disclination network for timescale of up to seconds. 
In physical units, the size of one unit cell is about 0.5 fim. 

In our simulations, a fraction (typically) of the simulated lattice was 
initialised in a BPII structure, previously equilibrated at the selected values 
of K, and r, whereas the rest of the lattice was initialised in the cholesteric 
phase. We first consider the case of low chirality and of a BPII domain 
growing inside a cholesteric phase (k = 1, Fig. [T]). It can be seen that the 
growth proceeds in a regular way: the BPII structure grows, as one would 
expect as its free energy is lower than that of the cholesteric (we are in a 
region in parameter space in which the cubic BPs are more advantageous 
structure than the cholesteric phase). The disclinations twist at the domain 
boundary to then merge into regular arrays of BPII unit cells. Their sizes 
are initially slightly larger than that of the initial template, but they rapidly 
equilibrate during growth to yield a virtually defect-free BPII lattice at the 
end of the simulations. 
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Increasing the chirality (Fig. [2]), we observe an interesting phenomenon. 
The growth is much faster, and is more irregular at the surface. The twist 
at the advancing boundary can no longer be accomodated within the BP II 
lattice. Therefore, for large enough values of the chirality, it becomes more 
advantageous to switch to a completely different blue phase lattice, with 
hexagonal structure. This hexagonally ordered structure is different from the 
previously reported hexagonal phases, which were theoretically predicted to 
be stable only in weak electric fields [35]. The BP we report here is not, 
so far as we know, observed experimentally, but it would be interesting to 
understand whether it could be related to local ordering in blue phase III, 
which appear at a larger chirality than either BPI or BPII. It is quite possible 
anyway for such a phase to arise as a metastable intermediate. 

The domain wall growth is very anisotropic. The results we have pre- 
sented are with the cholesteric helix pointing along the y axis, perpendicular 
to the domain wall. We have performed additional simulations with different 
relative orientation between helix and domain wall. If the helix lies on the 
plane of the domain wall, we do not observe the transition to the hexago- 
nally ordered structure. The anisotropy is also clear in the orientation we 
have chosen (Figs. [T]and|2]), as the speed of advance of the domain boundary 
is not equal at both ends. 

We have also performed simulations in which BPII was growing inside an 
isotropic phase at the same values of the reduced temperature. Again we find 
evidence of a transition between different BP lattices. At high values of the 
chiralities the "hexagonal" BP again forms, albeit through a different kinetic 
mechanism: instead of forming twisted disclination loops which are joined 
later on as within the cholesteric phase, the structure extends more regularly 
(Fig. |3]). This reinforces the suggestion that at high chirality the hexagonal 
ordering of the BP lattice should be thermodynamically advantageous. 

Finally, while Figs. [T} [2| [3] show the time evolution of the disclination 
line network. Fig. [i] gives the corresponding predicted optical patterns (one 
should keep in mind that the achievable resolution in an experiment would 
be smaller). It can be seen that the hexagonally ordered phase leaves a 
distinct signature from the BPII structure. Note that the hexagonal ordering 
is particularly evident in the larger corner spots in the optical pattern as well 
as generally in the -\/3 : 1 aspect ratio of the rectangular pattern. 
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4. Conclusions 

In conclusion, we have presented a hybrid lattice Boltzmann algorithm to 
solve the Beris-Edwards hydrodynamic equations of motion of cholesteric 
blue phases. The hybrid methodology mainly allows us to decrease the mem- 
ory requirements, which in previous full LB approaches eventually limited the 
lattice sizes which could be reached. 

We have shown that it is possible with this algorithm to perform suprado- 
main simulations of blue phases, in which a large number of unit cells is 
simultaneously considered, as opposed to the previously reported single unit 
cell simulations, which assumed a regular lattice of disclination lines. In par- 
ticular we have reported results on the growth of blue phase II domains in 
cholesterics, and shown that the growth dynamics is highly non-trivial. For 
low values of the chirality, the growth is regular and the blue phase forms 
without defects, whereas for large values of the chirality the growing region 
no longer has the structure of blue phase II, but attains a hexagonal ordering, 
and several defects appear at late times when the blue phase structure has 
invaded the whole computational domain. Finally, for very large values of 
the chirality, we observe a more regular hexagonal blue phases at late times. 
The switch between blue phase II and this hexagonal structure is due, in our 
interpretation, to the increasing twist which is developed at the domain wall, 
and which can no longer relax into a regular blue phase II unit cell. Our 
predicted transition from blue phase II to hexagonal structure appears for 
rather large values of the chirality, and it would be interesting to see whether 
these could be reached by the new temperature-stabilised blue phases. Also, 
it would be of interest to see whether this transition can be explained with 
extensions of the existing semi-analytical theories of blue phases, which are 
based on truncated expansions of the order parameter in Fourier components. 

We acknowledge EPSRC grant EP/E045316/1 for funding, and. G. P. 
Alexander and E. Orlandini for useful discussions and cpu time on Hector 
funded by EP/F054750/1. 
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Figure 1: (a)-(f) Snapshots of the 2D projection of the dischnation network of a growing 
BPII domain inside a cholesteric phase. Times (in simulation units) are shown in each of 
the panels. The chirality and reduced temperature were k — 1 and t — respectively. 
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Figure 2: (a)-(f) Snapshots of the 2D projection of the dischnation network of a growing 
BPII domain inside a cholesteric phase. Times (in simulation units) are shown in each of 
the panels. The chirality and reduced temperature were k — 2 and t — Q respectively. 
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Figure 3: (a)-(f) Snapshots of the 2D projection of the dischnation network of a growing 
BPII domain inside an isotropic phase. Times (in simulation units) are shown in each of 
the panels. The chirality and reduced temperature were k — 2 and t — respectively. 
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Figure 4: Predicted optical patters for the (a) initial configuration, and final simulated 
configurations for (b) k = 1 and a BPII in an initially cholesteric phase; (c) k = 2 and a 
BPII in an initially cholesteric phase; (d) k = 2 and a BPII in an initially isotropic 
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